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ABSTRACT 


An  elementary,  but  useful  problem  In  reliability  design  is  used  to 

Illustrate  some  of  the  principles  and  concerns  of  computational  probability. 

A  number  n  of  items  with  exponential  life  times  are  placed  in  parallel. 

The  time  until  all  items  have  failed  is  denoted  by  U  .  It  is  desirable  that 

n 

U  exceed  the  duration  Tf  where  T  is  a  random  variable.  Independent  of 

Q 

the  n  life  times.  We  wish  to  determine  the  smallest  value  of  n  for  which 
P(Un  >  T}  exceeds  1  -  e.  It  is  shown  that  if  T  has  a  delayed  distribution 
of  phase  type,  this  may  be  done  by  a  recursive  algorithm,  which  avoids  numer¬ 
ical  Integrations. 
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•  Formulation  of  the  Problem 

The  concern  for  algorithmically  tractable  solutions  to  probleas  In 
probability  adds  an  interesting  new  dimension  to  their  analysis.  In  the  con¬ 
struction  of  efficient  and  stable  algorithms,  certain  principles,  which  have 
not  been  given  much  attention  In  the  classical  analytic  presentations,  need 
to  be  followed.  In  this  paper,  some  of  these  principles  are  illustrated  by 
a  simple  example  from  reliability  theory. 

He  first  consider  n  items  with  independent,  identically  distributed, 
exponential  life  times  with  parsaeter  y.  These  items  are  placed  In  parallel 

and  system  failure  occurs  at  the  time  0  when  all  n  items  have  expired. 

u 

The  parallel  system  is  to  be  used  during  a  mission  of  duration  T,  where  T 
is  a  random  variable,  independent  of  the  life  times  of  the  items.  The  dis¬ 
tribution  of  T  is  denoted  by  F(«). 

The  probability  q(n)  -  P{UQ  <_  T>,  that  the  system  falls  before 
the  end  of  the  mission  is  readily  given  by 

m 

(1)  q(n)  -  |  (l-e’***)*1  dF(x)  ,  for  n>l. 

0 

Clearly  q(n)  is  strictly  decreasing  in  n  and  tends  to  saro  as  n  tends 
to  infinity.  For  every  e  >  0,  the  quantity  a*  ■  min  {n:  q(n)  <  c>  is 
well-defined,  n*  is  the  smallest  number  of  items,  which  will  yield  a  re¬ 
liability  of  at  least  1  -  c  . 

The  numerical  computation  of  n*  requires  some  care.  It  would  e.g. 
be  inadvisable  to  rewrite  (1)  as 

(2)  q(n)  -  <-l)k  (“]  f(kw>  , 

where  f (•)  is  the  Laplace-Stieltjes  transform  of  F(*)>  For  large  n,  the 
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quantity  q(n)  la  given  by  (2)  a a  tha  sun  of  tana  of  alternating  sign. 


Nuaarlcal  Imp  lament  Tina  of  (2)  uaually  auffara  from  loaa  of  significance. 
Ibis  servea  to  illustrate  our  first  point.  Mathematically  aqulvalant 


putt Ion,  k  number  of  tha  published  analytic  solutions  to  stochastic  nodals 
to  worthless  for  or  hasardous  In  nuaarlcal  laplaasntatlon.  Algorithmic 


tractablllty  deserves  to  be  Included  In  technical  discussions  of  probability 
aodals  to  a  far  greater  extent  than  is  presently  done. 

In  evaluating  the  quantities  q(n)  by  nuaarlcal  integration  In  Formula 
Cl),  considerable  cars  la  needed.  For  large  valuee  of  n,  the  Integrand  is 
nail  for  snail  to  aodarate  x.  The  laproper  liemann-StieltJes  Integrals  need 


to  be  truncated  at  a  high  value  of  x,  so  as  to  avoid  neglecting  a  substantial 
tall  contribution.  For  specific  probability  distributions  F(*) ,  we  nay  often 
perform  an  appropriate  change-of-varlable,  so  that  e.g.  the  Interval  of  inte¬ 
gration  will  bacon  finite  or  so  that  a  quadrature  method  may  be  applied.  This, 
of  course,  depends  on  particular  features  of  F(>). 

In  any  case  -  and  this  Is  our  second  point  -  It  will  not  be  possible  to 
write  one  computer  code  to  handle  a  very  wide  variety  of  probability  distri¬ 


butions  F(«).  Except  for  a  snail  nuaber  of  classical  particular  probability 
distributions,  such  as  e.g.  the  gasan  family,  it  Is  necessary  to  analyse  the 
numerical  Integration  procedure  from  first  principles.  This  is  not  a  matter 
of  mere  computer  programming.  For  this  -  and  a  fortiori  for  more  complex 
problems  -  this  analysis  should  be  the  concern  and  responsibility  of  the 


probab 111st. 


In  mny  problems,  we  may  exploit  particular  probabilistic  features  to 
expedite  the  computations.  For  the  versatile  class  of  distributions  of  phi 
type,  we  shall  show  that  the  quantities  q(n)  may  be  evaluated  without  nm 


Integrations .  Before  doing  so,  ve  shell  consider  a  more  general  version  of 
the  problen.  Let  there  be  n  iteas  originally,  which  fall  according  to  a 
pure  death  process,  i.e.  the  tines  between  successive  failures  are  independent 

and  exponentially  distributed  with  paraaeters  . The  tine -to - 

failure  then  has  a  generalised  Erlang  distribution.  If  (t)  denotes 
the  probability  that  k  of  die  n  original  iteas  are  still  alive  at  tins  t, 
the  quantity  q(n)  is  given  by 

m 

(3)  q(n)  -  J  ^n)(x)dF(x)  ,  for  n>0. 

0 

The  earlier  problea  corresponds  to  the  choice  ■  Jp,  for  1  <_ J  <_ n. 


2.  The  Case  Where  F(Q  is  of  Phase  Type 

F(»)  is  a  distribution  of  phase  type  (PH)  if  and  only  if  there  is  e 
finite-state  Markov  process,  with  one  absorbing  state  and  all  other  states 
transient,  in  which  F(*)  is  the  distribution  of  the  tlae  till  absorption 
[1,2].  The  generator  Q  of  such  a  Markov  process  is  then  of  the  fora 


Q  • 


T 

0 


I* 

0 
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where  T  is  a  square  aatrlx,  say  of  order  a,  with  negative  diagonal  ela- 
aents  and  nonnegative  off-diagonal  eleasnts.  The  aetrix  T  is  nonsingular 
and  T*  +  r*  ■  0  .  The  initial  probability  vector  is  denoted  by  [0,0^^] . 
The  initial  FH -distribution  F(>)  is  then  given  by 


(4) 


F(x)  •  1  -  o  exp  (Tx)  e  ,  for  *  >  0  . 


In  the  sequel,  we  shall  assuae  for  convenience  that 
F(>)  is  continuous  on  to,*). 


•  0,  so  that 


Let  now  ♦^(n)  be  eh*  conditional  probability  that  UQ  <_  T,  given  that 
th*  Markov  process  Q  start*  In  tha  state  1.  The  *- vector  with  components 
♦±(n) ,  1  <  i  <  s,  la  denoted  by  j£(n) .  It  Is  then  clear  that 

(5)  q(n)  »  a  i(n)  ,  for  n  >_  0  . 

the  vectors  £(j),  0  <_  j  <  n,  are  recursively  given  by 

(6)  itO)  -  *  , 

KJ)  -  VjCWjI  -  T)"1  l(J-l)  ,  for  1  <  j  <  n  . 

It  la  clear  that  £(0)  ■  *_.  Furthermore,  by  conditioning  on  the  tine  of  the 
first  of  the  j  failures,  given  that  th*  systan  starts  with  J  >_  1  itens, 
w*  obtain  by  using  tha  Markov  property  that 

7 

iCj)  -  J  exp  (Tx)  e  *  dx  •  i(j-l)  , 

0 

which  readily  yields  (6). 

From  Formulas  (5)  and  (6),  we  see  that  th*  q(n)  nay  be  ccnputed 
by  th*  numerical  solution  of  systens  of  linear  equations.  He  evaluate  th* 
successive  vectors  i(j),  J  >  1,  and  stop  as  soon  as  a£(J)  becomes  less 
then  th*  prescribed  c.  The  Integer  n*  Is  then  given  by  th*  Index  j  at 
which  computation  halts. 

Th*  numerical  solution  of  the  linear  equations 

(7)  i(J)  -  vj1  T  i(J)  -  i(J-l)  , 

is  highly  stable.  In  particular,  they  are  Ideally  suited  for  Gauss-Seldel 
Iteration.  As  we  shall  discuss  below,  there  nay  be  some  merit  In  using  this 
procedure,  although  It  Is  somewhat  slower  then  direct  Gauss  elimination. 
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3.  The  Case  Where  F(Q  is  a  Delayed  PH-dlstrlbutlon. 

In  practice,  it  la  unrealistic  to  assume  that  the  distribution  F(*) 
has  support  on  (0,“).  The  duration  T  will  usually  exceed  some  fixed  value 
a  >  0  with  probability  one.  It  Is  easy  to  modify  the  preceding  algorithm  to 
handle  this  case  also.  Let  F(*)  be  of  the  fora 


F(x)  *  0  ,  for  x  <,  a  , 


“l-o  exp  [T(x-a>]  e  ,  for  x^a  . 


It  may  be  called  a  delayed  PH-diatribution. 

The  quantity  q(n)  Is  then  given 

(8)  q(n)  -  1  ♦5“^ (a)  4(k)  ,  for  n>0  , 

k-0 

where  the  4(k)  are  the  quantities,  given  by  Formula  (5)  .  In  the  particular 
model,  with  which  we  started  this  discussion.  Formula  (8)  becomes 


q(n)  •  Z  [?]  (l-e-*1*)  e  4(k)  ,  for  n  >  0  . 

k-0  W 

The  binomial  probabilities  ^n^(a)  are  easily  evaluated.  In  general,  the 
probabilities  ♦^(t)  may  be  computed  by  the  numerical  solution  of  the  simple 
differential  equations 


(9) 

3t  *ka)<t)  ■  •*,  +  “irt-i  • 


for  1  <  k  <  n  , 
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for  t  >  0,  with  the  initial  conditions  $^(0)  *  1,  $5n^(0)  -  0,  for 

n  k 

0  <_  k  <  n. 

We  believe  that  the  preceding  discussion  illustrates  the  point  that  for 
many  probleas  of  practical  interest,  probabilistic  insight  nay  lead  to  natural 
and  efficacious  algorithms.  There  is  one  additional  point  to  be  made  in 
conclusion.  Our  simple  example  again  serves  to  illustrate  that  point  well. 

In  practical  situations,  we  often  wish  to  vary  parameters  of  the  model 
in  a  systematic  manner.  It  is  desirable  to  "nest"  the  successive  imolementa- 


If,  for  example,  we  wish  to  vary  the  parameter  a  in  Formula  (8),  it 
suffices  to  start  with  the  smallest  value  of  a  and  to  Integrate  the  dif¬ 
ferential  equations  (9)  progressively  up  to  the  largest  desired  value  of  a. 

The  sane  sequence  400 ,  k  >  0,  is  used  throughout. 

If  we  wish  to  modify  one  or  more  elements  of  the  matrix  T,  we  can  solve 
the  equations  (7)  by  an  iterative  method  and  use  the  preceding  solution 
vectors  as  starting  solutions. 

These  considerations  are  particularly  important  when  the  actual  computa¬ 
tions  are  carried  out  in  die  conversational  mode  et  a  remote  terminal.  The 
algorithm  then  becomes  a  powerful  tool,  which  may  be  "interrogated"  to 
assess  the  effect  of  varying  the  parameters  of  the  model.  In  complex  models, 
it  is  mostly  impossible  to  do  this  by  analytic  methods  only.  The  tines  between 
the  successive  printouts  of  results  sre  usually  a  dead  loss  to  the  user  and 
it  is  desirable  to  keep  these  as  short  as  possible. 

Once  the  algorithmic  analysis  of  one  stochastic  model  has  been  carried 
out,  it  often  provides  insights  into  the  analysis  of  related  models.  The 
preceding  discussion  may  easily  be  modified  to  handle  the  case  where  system 
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failure  occurs  when  only  r  out  of  the  original  n  items  survive.  The 
times  between  successive  failures  may  also  be  chosen  to  have  distributions 
of  phase  type  Instead  of  the  (particular)  exponential  distributions ,  used  In 
our  discussion.  The  algorithm  then  becomes  more  complex,  but  the  underlying 
mathematical  Ideas  remain  the  same.  To  Identify  the  structural  properties, 
useful  In  the  construction  of  efficient  algorithms.  Is  the  proper  objective 
of  computational  probability. 
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